!// Dmitry.Konovalov@jcu.edu.au Dmitry.A.Konovalov@gmail.com, 2013

module grid_mod; 
use memory_mod;
implicit none;
logical, save :: GRID_TRACE = .TRUE.;

private;
! NOTE 9-Apr-2013: name any new public functions as grid_xxx()
public test_grid_mod &
  , grid_lcr, grid_lcr_xstep, grid_lcr_wf, grid_lcr_init, grid_lcr_get_pow_r, GRID_LCR_T &
  , grid_step, grid_step_assert_4 ;


! grid transform: f(x) = r, where r>=0
!// x = ln(CR), CR=c+r; c > 0
!// exp(x) = c + r;  r = exp(x) - c
type GRID_LCR_T
  real(DP) :: rmin, rmax;
  real(DP) :: xmin, xmax, c, xstep;
  real(DP), pointer :: rv(:), r2(:), xv(:);  ! r, r^2, x
  real(DP), pointer :: cr(:), cr2(:);  ! (c+r), (c+r)**2
  real(DP), pointer :: crDivR(:), divR(:), sqcr1(:);  ! cr/r, (cr/r)**2, 1/r, 1./sqrt(c+r)
  real(DP), pointer :: xqdr(:);! x-quadr's
  real(DP), pointer :: xqdr_cr(:), xqdr_cr2(:);! xqdr_cr = xqdr*cr; xqdr_cr2 = xqdr*cr*cr;
  real(DP), pointer :: xqdr_cr2DivR(:), xqdr_cr2DivR2(:);! xqdr*cr**2/r, xqdr*(cr/r)**2;
  integer :: nx;
  
!  h == xstep
!  h2 = h / 2.0;
!  h3 = h / 3.0;
!  h90 = h / 90.;
!  h114 = 114.0 * h90;
!  a34 = 34. * h90;
  real(DP) :: h2, h3, h90, h114, a34;
end type


contains 

SUBROUTINE test_grid_mod();  
  character(LEN=64) :: name = "test_grid_mod";
  call dbg_in(name);
  call test_step_grid();
  !ok = vec_dot_test_speed();
  call dbg_out(name); 
end SUBROUTINE


SUBROUTINE test_step_grid();  
  integer :: nv;
  real(DP) :: minv, maxv, vstep;  
  real(DP), pointer :: pv(:), sg(:); ! sg-step grid
  character(LEN=64) :: name = "test_step_grid";
  call dbg_in(name);

  nv = 2;  minv = D0; maxv = D1;
  pv => mem_vec(nv);  
  sg => mem_vec(nv);  
  call grid_step(pv, minv, maxv, vstep);
  call dbg_assert_zero("2.pv(1)=minv", pv(1) - minv);
  call dbg_assert_zero("2.pv(nv)=maxv", pv(nv) - maxv);

  call grid_from_step(sg, minv, vstep);
  call dbg_assert_zero("pv(1)==sg(1)", pv(1) - sg(1));
  call dbg_assert_zero("pv(2)==sg(2)", pv(2) - sg(2));
  call mem_del_vec(pv);
  call mem_del_vec(sg);

  nv = 3;  minv = -D1; maxv = D1;
  pv => mem_vec(nv);  
  sg => mem_vec(nv);  
  call grid_step(pv, minv, maxv, vstep);
  call dbg_assert_zero("3.pv(1)=minv", pv(1) - minv);
  call dbg_assert_zero("3.pv(2)=0", pv(2));
  call dbg_assert_zero("3.pv(nv)=maxv", pv(nv) - maxv);

  call grid_from_step(sg, minv, vstep);
  call dbg_assert_zero("pv(1)==sg(1)", pv(1) - sg(1));
  call dbg_assert_zero("pv(2)==sg(2)", pv(2) - sg(2));
  call dbg_assert_zero("pv(3)==sg(3)", pv(3) - sg(3));

  call mem_del_vec(pv);
  call mem_del_vec(pv); ! testing; call for a freed ptr
end SUBROUTINE

SUBROUTINE grid_step(v, vmin, vmax, step_out);      
  real(DP), intent(OUT) :: v(:);      
  real(DP), intent(OUT), optional :: step_out;
  real(DP), intent(in) :: vmin, vmax; ! v1-first value; v2-last value
  integer :: mini, maxi, i, nv; ! nv - num of point in v
  real(DP) step;   
  call dbg_str("-->grid_step()", GRID_TRACE);
  call vec_range(mini, maxi, v);  
  nv = maxi - mini + 1;
  call dbg_assert_true("(nv > 0)=", nv > 1);
  step = (vmax - vmin) / real(nv-1, DP);
  if (present(step_out))  step_out = step;
  v(mini) = vmin;
  do i = mini + 1, maxi
    v(i) = vmin + step * real(i-mini, DP);
  end do
end SUBROUTINE

SUBROUTINE grid_from_step(v, vmin, vstep, vmax_out);      
  real(DP), intent(OUT) :: v(:);      
  real(DP), intent(OUT), optional :: vmax_out;
  real(DP), intent(in) :: vmin, vstep; ! v1-first value; v2-last value
  integer :: mini, maxi, i, nv; ! nv - num of point in v
  !call dbg_str("-->grid_step()", GRID_TRACE);
  call vec_range(mini, maxi, v);  
  nv = maxi - mini + 1;
  call dbg_assert_true("(nv > 0)=", nv > 1);
  v(mini) = vmin;
  do i = mini + 1, maxi
    v(i) = vmin + vstep * real(i-mini, DP);
  end do
  if (present(vmax_out))  vmax_out = v(maxi);
end SUBROUTINE


! only checks two intervals at the front and two at the back
SUBROUTINE grid_step_assert_4(px);     
  real(DP), intent(IN) :: px(:);      
  real(DP) st, st2;
  integer xn, mini, maxi;
  !call dbg_str("-->grid_step_assert_4()", GRID_TRACE);
  call vec_range(mini, maxi, px, xn);
  call dbg_assert_true("xn >= 3", xn >= 3); ! pointless if < 3
  call dbg_assert_true("px(mini+1)>px(mini)", px(mini+1) > px(mini));
  call dbg_assert_true("px(maxi)>px(maxi-1)", px(maxi) > px(maxi-1));
  
  st = px(mini + 1) - px(mini);
  st2 = px(mini + 2) - px(mini + 1);
  call dbg_assert_abs_diff_LE("mini", st, st2, EPS14);
  st = px(maxi) - px(maxi - 1);
  call dbg_assert_abs_diff_LE("mini,maxi", st, st2, EPS14);
  st2 = px(maxi - 1) - px(maxi - 2);
  call dbg_assert_abs_diff_LE("maxi", st, st2, EPS14);
end SUBROUTINE


! LCR GRID

!  // f(x) = ln(CR), CR=c+r; c > 0
SUBROUTINE lcr_c_from_x0_r0(lcr);  
  type(GRID_LCR_T), intent(INOUT) :: lcr;
  call dbg_dp("lcr%xmin", lcr%xmin, GRID_TRACE);
  call dbg_dp("lcr%rmin", lcr%rmin, GRID_TRACE);
  lcr%c = EXP(lcr%xmin) - lcr%rmin;
  call dbg_dp("lcr%c", lcr%c, GRID_TRACE);
  call dbg_assert_true("lcr%c > 0", lcr%c > 0);
end SUBROUTINE

!  // f(x) = ln(CR), CR=c+r; c > 0
real(DP) function r_to_lcr(r, lcr) result(x);  
  type(GRID_LCR_T), intent(IN) :: lcr;    
  real(DP), intent(IN) :: r;
  x = LOG(lcr%c + r);
end function

!  // f(x) = r,
!  // x = ln(CR), CR=c+r; c > 0
!  // exp(x) = c + r;  r = exp(x) - c
real(DP) function lcr_to_r(x, lcr) result(r);  
  type(GRID_LCR_T), intent(IN) :: lcr;    
  real(DP), intent(IN) :: x;
  r = EXP(x) - lcr%c;
end function


SUBROUTINE vec_lcr_to_r(res, xv, lcr);  
  real(DP), intent(OUT) :: res(:);    
  real(DP), intent(IN)  :: xv(:);      
  type(GRID_LCR_T), intent(IN) :: lcr;
  integer :: i, mini, maxi;
  call dbg_str("-->vec_lcr_to_r()", GRID_TRACE);
  call vec_init_res(res, xv);
  call vec_range(mini, maxi, xv); 
  do i = mini, maxi 
    res(i) = lcr_to_r(xv(i), lcr);
  end do
end SUBROUTINE


SUBROUTINE grid_lcr_init(lcr);  
  type(GRID_LCR_T), intent(OUT) :: lcr;
  call dbg_str("-->grid_lcr_init()", GRID_TRACE);
  lcr%nx = 0;
  lcr%rmin = D0;  lcr%rmax = D0;
  lcr%xmin = D0;  lcr%xmax = D0;  lcr%c = D0; 
  lcr%xstep = D0; lcr%h2 = D0; lcr%h3 = D0; lcr%h90 = D0; lcr%h114 = D0; lcr%a34 = D0; 
  
  call mem_del_vec2(lcr%rv, lcr%r2);
  call mem_del_vec(lcr%xv);
  call mem_del_vec2(lcr%cr, lcr%cr2);
  call mem_del_vec(lcr%crDivR);
  call mem_del_vec(lcr%divR);
  call mem_del_vec(lcr%sqcr1);
  call mem_del_vec(lcr%xqdr);
  call mem_del_vec2(lcr%xqdr_cr, lcr%xqdr_cr2);
  call mem_del_vec2(lcr%xqdr_cr2DivR, lcr%xqdr_cr2DivR2);
end SUBROUTINE



! rwfv - radial wave function vector
subroutine grid_lcr_wf(fv, lcr);  
  real(DP), intent(INOUT) :: fv(:); ! fv - function vector
  type(GRID_LCR_T), intent(IN) :: lcr;
  call vec_assert_math(fv, lcr%sqcr1);
  fv = fv * lcr%sqcr1; ! 1/sqrt(c+r)
end subroutine


function grid_lcr_get_pow_r(lcr, k) result(res);  
  real(DP), pointer :: res(:); 
  type(GRID_LCR_T), intent(IN) :: lcr;
  integer, intent(IN) :: k;
  if (k == 1) then
    res => lcr%rv;
    return;
  elseif (k == 2) then
    res => lcr%r2;
    return;
  endif
  call dbg_assert_true("grid_lcr_get_pow_r(k not ready", .FALSE.); !/stop here
end function


!// input: rmin, rmax, xmin, nx
SUBROUTINE grid_lcr(lcr);  
  type(GRID_LCR_T), intent(INOUT) :: lcr;
  character(LEN=64) :: name = "grid_lcr";
  call dbg_in(name, GRID_TRACE);
  
  call dbg_assert_true("lcr%nx > 1", lcr%nx > 1);
  call dbg_dp("lcr%xmin", lcr%xmin, GRID_TRACE);
  call dbg_dp("lcr%rmin", lcr%rmin, GRID_TRACE);
  call dbg_dp("lcr%rmax", lcr%rmax, GRID_TRACE);
  call lcr_c_from_x0_r0(lcr);
  lcr%xmax = r_to_lcr(lcr%rmax, lcr);
  call dbg_dp("lcr%xmax", lcr%xmax, GRID_TRACE);
  !call dbg_assert_true("lcr%xmax > 0", lcr%xmax > 0);

  lcr%xv => mem_vec(lcr%nx);
  call grid_step(lcr%xv, lcr%xmin, lcr%xmax, lcr%xstep);
  call grid_lcr_from_xv(lcr);  
end SUBROUTINE

!// input: rmin, xmin, xstep, nx
SUBROUTINE grid_lcr_xstep(lcr);  
  type(GRID_LCR_T), intent(INOUT) :: lcr;
  character(LEN=64) :: name = "grid_lcr";
  call dbg_in(name, GRID_TRACE);
  
  call dbg_assert_true("lcr%nx > 1", lcr%nx > 1);
  call dbg_dp("lcr%xmin", lcr%xmin, GRID_TRACE);
  call dbg_dp("lcr%xstep", lcr%xstep, GRID_TRACE);
  call dbg_dp("lcr%rmin", lcr%rmin, GRID_TRACE);

  lcr%xv => mem_vec(lcr%nx);
  call grid_from_step(lcr%xv, lcr%xmin, lcr%xstep, lcr%xmax);
  call lcr_c_from_x0_r0(lcr);
  lcr%rmax = lcr_to_r(lcr%xmax, lcr);
  
  call grid_lcr_from_xv(lcr);  
end SUBROUTINE

SUBROUTINE grid_lcr_from_xv(lcr);  
  type(GRID_LCR_T), intent(INOUT) :: lcr;
  real(DP) :: h;
  character(LEN=64) :: name = "grid_lcr_setup";
  call dbg_in(name, GRID_TRACE);

  lcr%rv => mem_vec(lcr%nx);
  lcr%r2 => mem_vec(lcr%nx);
  lcr%cr => mem_vec(lcr%nx);
  lcr%cr2 => mem_vec(lcr%nx);
  lcr%crDivR => mem_vec(lcr%nx);
  lcr%divR => mem_vec(lcr%nx);
  lcr%sqcr1 => mem_vec(lcr%nx);

  call vec_lcr_to_r(lcr%rv, lcr%xv, lcr);
  lcr%cr = lcr%rv + lcr%c; 
  lcr%r2 = lcr%rv * lcr%rv; 
  lcr%cr2 = lcr%cr * lcr%cr; 
  call vec_div_SAFE(lcr%divR, D1, lcr%rv);
  lcr%crDivR = lcr%cr * lcr%divR;
  lcr%sqcr1 = D1 / SQRT(lcr%cr);
  
  ! the following is needed in yk
  h = lcr%xstep;
  lcr%h2 = h / D2;
  lcr%h3 = h / D3;
  lcr%h90 = h / D90;
  lcr%h114 = D114 * lcr%h90;
  lcr%a34 = D34 * lcr%h90;
  
end SUBROUTINE

end module 
